High-sensitivity dynamic diffuse fluorescence tomography system for fluorescence pharmacokinetics

Abstract. Significance: Dynamic diffuse fluorescence tomography (DFT) can recover the static distribution of fluorophores and track dynamic temporal events related to physiological and disease progression. Dynamic imaging indocyanine green (ICG) approved by the food and drug administration is still under-exploited because of its characteristics of low quantum yield and relatively rapid tissue metabolism. Aim: In order to acquire the ICG tomographic image sequences for pharmacokinetic analysis, a dynamic DFT system was proposed. Approach: A fiber-based dynamic DFT system adopts square-wave modulation lock-in photon-counting scheme and series-parallel measurement mode, which possesses high sensitivity, large dynamic range, high anti-ambient light ability in common knowledge, as well as good cost performance. In order to investigate the effectiveness of the proposed system, the measurement stability and the anti-crosstalk—a crucial factor affecting the system parallelization—were assessed firstly, then a series of static phantoms, dynamic phantoms and in vivo mice experiments were conducted to verify the imaging capability. Results: The system has the limited dynamic range of 100 dB, the fluctuation of photon counting within 3%, and channel-to-channel crosstalk ratio better than 1.35. Under the condition of a sufficient signal-to-noise ratio, a complete measurement time for one frame image was 10.08 s. The experimental results of static phantoms with a single target and three targets showed that this system can accurately obtain the positions, sizes, and shapes of the targets and the reconstructed images exhibited a high quantitativeness. Further, the self-designed dynamic phantom experiments demonstrated the capability of the system to capture fast changing fluorescence signals. Finally, the in vivo experiments validated the practical capability of the system to effectively track the ICG metabolism in living mice. Conclusions: These results demonstrate that our proposed system can be utilized for assessing ICG pharmacokinetics, which may provide a valuable tool for tumor detection, drug assessment, and liver function evaluation.


Introduction
lock-in photon-counting technique for realizing simultaneous irradiation of multiple light sources and parallel detection of multiple detectors. To validate the effectiveness of the proposed system, a series of static and dynamic phantom experiments, as well as in vivo experiments were conducted.

Dynamic DFT System
As illustrated in Fig. 1, the proposed dynamic DFT system consists of a light source unit for generating continuous-wave laser light, a detection unit for detecting the excitation and emission lights, and a user interface unit for controlling and visualizing imaging results.
In this system, 16 pig-tailed laser diodes (L785P25, Thorlabs) at 785 nm wavelength driven by a constant current source module are controlled by square wave generator implemented on a field-programmable gate-array (FPGA, DE0-Nano-Soc Kit with Cyclone V, Terasic) hardware with the frequency modulation ranges from 0 to 10 kHz. Sixteen source fibers having a core diameter of 62.5 μm and a numerical aperture (NA) of 0.275 coupled with 16 detection fibers having a core diameter of 500 μm and a NA of 0.370 are mounted into the fiber holders with a uniform distribution around a cylindrical imaging chamber fixed on a motorized translation stage. The source fibers conduct the lights to illuminate the object fixed in the imaging chamber. Simultaneously the transmitted lights on the object boundary are collected by detector fibers and directed into an integrated module of programmable 4-1 × 4 fiber optic switch (FSW4-1 × 4-MM-C, Guilin Institute of Optical Communications, China). The light outputs from the fiber optic switch are collimated using collimators (FC230FC-B,Thorlabs) for normal incidence to four six-hole motorized filter wheels (FF01-832/37-25, Semrock) where one of the holes houses a single-band pass filter (FEL0700, Thorlabs) with the transmittance >95% at 830 AE 5 nm and <0.01% at 780 AE 5 nm for blocking excitation light. In principle, the fluorescence emission and excitation signals are detected with and without the filter stack in place, respectively. Subsequently, the light signals are fed into four PMT counting heads (H8259-02e, Hamamatsu, Japan) having a count linearity of 2.0 × 10 6 s −1 and a dark count of 400 s −1 , and thereby a limited dynamic range of 100 dB, and are converted into electrical pulses corresponding to the  PMT single-electron responses. Simultaneously, the collected composite signals are demodulated by digital lock-in photon-counting detection module. The whole setup is controlled by a LabVIEW-programmed graphical user interface. In this dynamic DFT system, lock-in photon-counting technique were employed to obtain high sensitivity, large dynamic range, and high ability to reject ambient light as a common knowledge. 21 In addition, serial-parallel detection mode was proposed to obtain a balance between measurement speed and cost-effectiveness. In the following, the detailed design principle was introduced.

Serial-parallel measurement mode
The light excitation and detection processes adopting serial-parallel measurement mode are shown schematically in Fig. 2. The detailed process for a complete measurement can be described as follows. First, four modulated light sources (nos. 1 to 4) illuminate the object simultaneously, and the other positions divided into three groups (nos. 5 to 8, nos. 9 to 12, and nos. 13 to 16) detect the transmission light in turn, by a 4-1 × 4 optical switch. In the same manner, nos. 5 to 8, nos. 9 to 12, and nos. 13 to 16 light sources are modulated, respectively, and the other three groups are taken as detection positions accordingly. Consequently, a complete measurement can acquire a total number of 16 × 12 detection data for one frame image reconstruction, and the measurement time can be calculated with the following equation E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 3 5 2 where T i is the integration time of digital lock-in photon-counting module; T s is the switching time of optical switch; the constants 3 and 4 denote the switch times of detection and source fibers, respectively.

Square-wave modulation lock-in photon-counting scheme
For enabling simultaneous multichannel excitation and detection, phase-lock-in technique with square wave modulation was employed and primarily implemented by the digital-phasedetection (DPD) blocks in the system. For each measurement, four laser diodes in a group are switched "ON" and "OFF" by the four square-wave modulation signals with a duty cycle of 50% at the frequencies of 3.2, 3.6, 4.0, and 4.4 kHz, respectively. Meanwhile, they are sent to the lockin photon-counting module as the reference signals of "1" weights for demodulation. Compared with sine-wave modulation scheme with the nonlinearity caused by the error of phase modulation, the square-wave modulation would not introduce nonlinearity error ascribe to only digital operations involved in the demodulation process, and more importantly, would greatly simplified the DPD block hardware. As illustrated in Fig. 3, each PMT-connected demodulator (with regards to a detection position) in the lock-in photon-counting module is configured with four DPD blocks to maximally discriminate the component signals at four modulation frequencies respectively, from the received composite light signal. The intensity of the demodulated component signal is achieved by the multi-periodic reference-weighted counting (RWC) strategy in the DPD block, as shown on the right of Fig. 3, where the reference weights at the time of the PMT outputs are accumulated for multiple periods. To eliminate the initial time-shift uncertainty of the signal, a pair of in-phase I ðsÞ ðtÞ and quadrature Q ðsÞ ðtÞ reference signals having the same frequency f ðsÞ (s ¼ 1, 2, 3, and 4) are introduced into each digital DPD block. These two reference signals are mathematically regarded as bipolar square-wave of "1" weights in the RWC operations, and can be expressed in terms of the Fourier series, which are given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 4 3 2 (2) The RWC in the digital DPD is mathematically equivalent to the multiplication-integration operations in the analog DPD, where the operations can be conducted by adding or subtracting "1" (multiplication) for multiple periods (integration) depending on the reference phase of the PMT pulse occurrences P ðsÞ d ðtÞ, as illustrated in Fig. 4 Data_ f (1) Reference_f (1) Output Block2 PMT_output Intergration _signal Output Reference_f (2) Output Reference_f (3) Reference_f (4) Data_ f (2) Data_ f (3) Fig. 3 Block diagram of digital-phase-detection based on reference-weighted counting strategy.
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 6 3 5 where A Before conducting imaging experiments, we have assessed the system performances in terms of stability and anti-crosstalk-a crucial factor affecting the system parallelization. The results showed that the photon counts fluctuated within 3% during 4 h period, and the channel-tochannel crosstalk ratio was better than 1.35%, suggesting the system demonstrated good performance. The detailed evaluation methods can be referenced in our previous work. 22 For the following experiments, the integration time of lock-in photon-counting module was determined to be 0.8 s to obtain a sufficient signal-to-noise ratio, and the switching time of optical switch was set to 0.04 s. Thus, a complete measurement time for one frame image was 10.08 s, according to Eq. (1).

DFT Reconstruction Method
NIR light transport in biological tissues can be described by the diffusion equation. For fluorescence tomographic imaging with continuous wave light source, the photon transport model can be described by a coupled diffusion equation, which is given by 23 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 3 2 8 where subscripts x and m denote the excitation and emission wavelengths, respectively; where I x ðr d ; r s Þ is the calculated excitation flux at r d for a source at r s ; Gðr d ; rÞ is the Green's function describing light propagation at the fluorescence wavelength. Equation (7) can be degenerated into a matrix equation I nb ¼ WX, where W is the weighted matrix of size N d×s × N voxels , where N d×s is the number of measurements and N voxels is the number of nodes of the finite element method meshes; X is the fluorescence yield to be reconstructed. In this paper, the fluorescence yield images were recovered relayed on algebraic reconstruction technique, where the iteration number and relaxation parameter were empirically determined to be 30 and 0.2, respectively. 23

Pharmacokinetic Analysis Method
Based on dynamic DFT measurement and reconstruction algorithm, a sequence of fluorescence yield images can be acquired. Further, the ICG concentration CðtÞ (μM) at time t can be calculated from fluorescence yield (ημ af ðtÞ) reconstruction by using the equation 24 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 6 1 4 μ af ðtÞ ¼ ln 10ξCðtÞ; (8) where ξ is the extinction coefficient of ICG (0.013 mm −1 μM −1 ). 25 For the in vivo experiment involving assessing liver function, the ICG kinetics in the livers can be monitored by measuring the ICG concentration-time curves that can be described by a biexponential fitting that is expressed as 26 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 5 3 3 CðtÞ ¼ −A expð−αtÞ þ B expð−βtÞ; (9) where A and B (unit: a.u.) are the zero-time intercepts, representing the initial hepatic ICG concentration; α and β (unit: min −1 ) are the uptake and excretion rates representing ICG influx to the liver and disappearance from the liver, respectively, which are used to quantitatively assess liver function. Figure 5 shows the sketches and photos of phantoms with a single target and three targets for static experiments. The phantoms were made of polyformaldehyde, which had the scattering characteristics similar to biological tissue. The absorption and reduced scattering properties were determined to be 0.0034 and 0.08 mm −1 at 785 nm, using time-resolved spectroscopy method. The diameter and height of the phantoms were 30 and 50 mm, respectively. To mimic a target tumor in biological tissue, a mixture of 1% Intralipid and ICG (Liaoning Tianyi Biological Pharmaceutical Co., Ltd., China) was placed in a circular hole drilled in the phantom. Here, two types of phantoms were designed: a single-target phantom with a hole of 5 mm diameter, 30 mm depth, and 7.5 mm off-center; a three-target phantom with three holes of 4 mm diameter and 35 mm depth, locating at the positions that were 9 mm away from the central axis of the phantom and 120 deg separating angles among them. In the single target experiment, the hole was filled with the solution of 1% Intralipid mixed with 0.5, 1.0, and 2.0 μmol∕L ICG, respectively. Furthermore, in order to verify the capability of imaging multiple targets simultaneously, two types of three-target phantoms were configured. One phantom was filled with 1% Intralipid mixed with 3.0 μmol∕L ICG in all three holes. The other was filled with 1% Intralipid mixed   Figure 6 shows a self-designed dynamic phantom. It was made of a cylindrical polyformaldehyde chamber of 30-mm diameter, filled with 1% Intralipid as the background and a cylindrical tube of 12-mm diameter filled with a mixture of 1% Intralipid solution and ICG as the target. To simulate the ICG concentration change process in the target region, 1% Intralipid fresh solution was pumped into the target tube, using a peristaltic pump-controlled silicone catheter (BT100-02, China). For evenly mixing the solution in the target cylinder during a limited time, the inlet catheter was placed at the height of 30 mm below the solution level (42 mm), and the imaging plane was selected at the height of 25 mm below but close to the inlet catheter.

Dynamic phantom experiments
According to the instruction of the peristaltic pump, the flow speed V pump (unit: mL/min) of the silicone catheter and the concentration of ICG can be calculated using the equation E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 3 1 7 where V pump (rpm) denotes the rotation speed of the peristaltic pump and α is a constant determined by the silicone catheter type. C 0 and V 0 are the initial concentration and volume of ICG, respectively. Under the same V pump , the change rate of ICG concentration can be altered by setting different initial ICG concentration. For performing dynamic experiment, 1% Intralipid mixed with 1.0 μmol∕L ICG was adopted as the initial target concentration. The rotation speed of peristaltic pump was adjusted to 4 rpm, and α was select to be 0.00105, which means the injection speed of the solution was 0.0252 ml∕ min. The whole detection process lasted ∼18 min for 104 complete measurements.  hair by using depilatory cream for better detection of light signal. Then a bolus of ICG (0.6 to 0.7 mL, 50 μg∕mL) was injected via tail vein. Subsequently, the mouse was fixed in the imaging chamber and filled the interspace with 1% Intralipid. According to anatomical principle, we ascertained the imaging plane of mouse liver at 0.5 to 1 cm blow the xiphoid and slightly marked the corresponding chamber surface. It is worth mentioning that the uptake rate of ICG in mice is rapid enough, which usually takes about 3 min to accumulate to the peak of ICG content in liver after injection. Thus, the dynamic probe started 3 min later after the ICG injection. In order to capture the metabolic process of ICG in liver, 120 complete measurements were implemented continuously, ∼20 min. To verify the reliability of liver position obtained by DFT reconstruction, a micro-CT imaging system (SkyScan 1276, bruker, German) was used to obtain the anatomical structure. Before micro-CT imaging, an iron anchor point was plastered on the previous mark position of the chamber surface to guarantee the same imaging plane for DFT and micro-CT. The micro-CT imaging began 30 min after intravenous administration of a contrast agent at a dose of 0.4 mL∕mouse, produced at Tianjin Medical University. It is worth mentioning that the bottom of the imaging chamber is detachable. Thus, the tail vein injections of ICG and liver contrast agent were able to be performed in the same chamber. For in vivo imaging, the μ av ðrÞ and μ 0 sv ðrÞ are set to 0.03 and 1.0 mm −1 , respectively, corresponding with the average values of bulk tissue. 27 3 Results and Discussions Figure 8(a) illustrates the reconstructed yield images in imaging plane and the corresponding x-profiles at y ¼ 0 mm. We can observe that the reconstructed target was clearly distinguishable from the background and displayed a good agreement with the true target (real line) in terms of the localization and size. We can also see that the reconstructed target region became more remarkable with the ICG concentration increase, and the x-profile plotted along the horizontal line passing through the centers of the target quantitively depicted that the fluorescence yield increased in proportion to ICG concentration. Figure 8(b) shows the reconstructed yield images, and the profiles along the lines passing through the centers of the phantoms and the targets. From the reconstructed fluorescence yield images of the three-target phantoms and the corresponding target profiles, we can observe that the proposed system can still disclose the locations and sizes of the three targets for both types of phantoms. For the phantom with three targets having the same ICG concentration, the reconstructed fluorescence yields of three targets demonstrated a good consistency in quantitativeness, except that the target 3 is slightly larger. For the three targets with different ICG concentration, it is evident that the reconstructed fluorescence yield increases with the increase of ICG concentration. In addition, we find that the fluorescence yield ratio calculated using the maximum of the target-profile approximately conforms to 1:2:4, which is consistent with the theoretical value. These results demonstrate that the system possesses high sensitivity and quantitativeness. We can see that the reconstructed target is consistent with the true location and size of the target (black circle). Figure 9(b) shows the theory curve (real line) of the ICG concentrations, the raw data (markers) calculated by the average ICG concentrations in the target region and the fitted curve (dotted line) versus the measurement time. It can be seen that the reconstructed ICG concentration demonstrated a declining trend with time, and the reconstructed values became closer to theoretical values with the measurement time increased, indicating the capacity of the dynamic DFT system to capture the rapidly changing fluorescence signals. However, the reconstructed values partly deviated from the theory values, which is mainly because that the exchange of the solution in the target needed some time. The fluorescence yield values in the dynamic phantom experiments could not decrease as fast as the theoretical calculating values, especially for the situation that the ICG concentration is higher.

In Vivo Experimental Results
As mentioned previously, the dynamic measurement started after 3 min post-injection. Thus, only the pharmacokinetic parameters of B and β were extracted based on the second exponential part of Eq. (9). Figure 10(a) shows a slice of the micro-CT images corresponding to the DFT measurement plane and typically illustrates the overlays of the anatomical images and the corresponding fluorescence tomographic image sequences with an interval of 14 image frames. For reflecting the average metabolic process of ICG in livers, the mean values of the ICG concentrations in liver were calculated. Figure 10(b) shows the plots of time-course changes in the normalized ICG concentrations and the corresponding fitted curves that can be obtained by an exponential-curve-fitting method. As can been seen from Fig. 10(b), the ICG concentrations decrease with the ICG metabolism in liver and the fitted curves are close to each other. To quantitatively assess the pharmacokinetic parameters of each mouse, the fitted parameters (B and β) for each mouse are listed in Table 1. As shown in Table 1, all mice show close values of Bð0.9758 AE 0.03513Þ and βð0.04558 AE 0.005207 min −1 Þ, suggesting that the pharmacokinetic parameters of ICG in liver exhibit a remarkable consistency within a specific kind of mice. The slight difference of the pharmacokinetic parameters may be caused by different experimental conditions and physiological status of the mice.

Conclusions
We proposed a dynamic DFT system for achieving pharmacokinetic parameters of ICG in living mice. Under the consideration of the balance between cost and effectiveness, the system employs lock-in photon-counting detection of square wave modulation mode to provide high sensitivity, large dynamic-range, and high ability to reject ambient light. The comprehensive assessments of the system have been experimentally conducted on static phantoms, dynamic phantom, and in vivo mice. The static phantom with a single-target and three-target experiments validate the capability of the system to recover the locations and sizes of the targets, and demonstrate an  excellent quantitativeness of the reconstructed fluorescence yield. Furthermore, the self-designed dynamic phantom experiments verified that the dynamic system can capture rapidly changing fluorescence signals. Finally, the in vivo experiments show the feasibility to acquire the timecourse of ICG concentration images in health mice's livers for pharmacokinetic analysis. These results demonstrate that our proposed system may provide a valuable tool for tumor detection, drug assessment, and liver function evaluation.

Disclosures
The authors declare that there are no conflicts of interest related to this article.